CS448 - Lab 2: Filter Design and Usage¶

For this lab you will learn how to design some simple filters and how to apply them to solve some common audio problems. Python’s scipy.signal package has an extensive set of commands to help you design filters (firwin, firwin2, butter, cheby1, cheby2, ellip, …), so there is no shortage of options.

Part 1: When to use what¶

There will be four cases to this part, each requiring a different type of filter to address a problem in a recording. The input recordings are:

  • case1.wav [https://drive.google.com/uc?export=download&id=1eTsDiXqqLQv3murPz25O54E89i3DL__d ] : A noise-corrupted speech signal. We want to eliminate the noise.
  • case2.wav [https://drive.google.com/uc?export=download&id=1egd22CxPUe6sINIi0FPTbMfG4S_In2hT ] : Same as above, different type of noise. We want to remove the noise again.
  • case3.wav [https://drive.google.com/uc?export=download&id=1eF-VOVWoT1rh1wAC06WT1ANusMKyDYSn ] : Bird songs during a thunderstorm. As a world renowned ornithologist you need to have a cleaner recording of the bird songs for further analysis.
  • case4.wav [https://drive.google.com/uc?export=download&id=1eeizGhrBICf6pW5OXcbq7ChF4m2N6yIk ]: The signal that we require to extract here is a Morse code which is buried in environmental noise. Design a filter to bring out the beeps.

For each case do the following:

  • Plot the spectrogram of the given sound and identify the problem
  • Describe what kind of filter will address this problem and why
  • Design an FIR filter using scipy.signal.firwin() and/or scipy.signal.firwin2()
  • Design an FIR filter using the formulas in my slides (i.e. do not use scipy.signal)
  • Design an IIR Butterworth filter using scipy.signal's routines
  • Show a plot comparing the response of all the filters (hint: scipy.signal.freqz)
  • Plot spectrograms of the filtered recordings and compare with the input to show that it fixed the problem
  • Play the filtered sounds so that we can hear how well it works

Make some observations on how the results differ between an FIR and IIR filter and try to find the best possible filter size/type/parameters to produce the best result. Show results under various parameters (e.g. filter length) and make some plots that demonstrate the effects of these parameters. Most importantly, try to get a sense of how these design choices impact audible quality. Being able to listen at a sound and identify what’s wrong and how to fix it is a big part of audio processing.

Hint: To apply an FIR filter you can use scipy.signal.convolve, to apply an IIR filter (or an FIR) you can use scipy,signal.lfilter.

In [ ]:
import matplotlib.pyplot as plt
import numpy as np
from scipy.io import wavfile
from scipy.signal import firwin, freqz, lfilter, get_window, butter, convolve
In [ ]:
# Make a sound player function that plays array "x" with a sample rate "rate", and labels it with "label"
def sound(x, rate=8000, label=''):
    from IPython.display import display, Audio, HTML
    display(
        HTML('<style> table, th, td {border: 0px; }</style> <table><tr><td>' +
             label + '</td><td>' + Audio(x, rate=rate)._repr_html_()[3:] +
             '</td></tr></table>'))

# Make a function that plots the spectrogram of a sound
def plot_spectrogram(input_sound, fs, title="Spectrogram"):
    plt.title(title)
    plt.specgram(input_sound, Fs=fs, cmap="winter")
    plt.xlabel('Time [sec]')
    plt.ylabel('Frequency [Hz]')
    plt.show()

# Make a function that plots the frequency response of a filter
def plot_freq_response(w,
                       h,
                       fs,
                       cutoff,
                       title):
    plt.plot(0.5 * fs * w / np.pi, np.abs(h), 'b')
    plt.axvline(cutoff, color='k', linestyle='--')
    plt.xlim(0, 0.5 * fs)
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.title(title)
    plt.grid(True)
    plt.show()
In [ ]:
# Load the sound, play it, and plot it's spectrogram
fs_case1, input_sound_case1 = wavfile.read("./data/case1.wav")
sound(input_sound_case1, rate=fs_case1, label='case1.wav')

plot_spectrogram(input_sound_case1, fs_case1, title="Spectrogram of case1.wav")

# Design a filter that fixes things
### A low-pass filter with a cutoff frequency around 6KHz would be good here as the
### noise present in the recording seems to be the only data above that threshold. I chose
### the exact frequency 5700Hz after a process of trial and error.
# Filter using firwin and filtering the input sound
firwin_cutoff_case1 = 5700
firwin_numtaps = 5000

firwin_case1 = firwin(firwin_numtaps,
                      firwin_cutoff_case1,
                      fs=fs_case1,
                      pass_zero=True)
firwin_sound_case1 = convolve(firwin_case1, input_sound_case1)

# Plot the resulting spectrogram, frequency response, and play the resulting sound
sound(firwin_sound_case1,
      rate=fs_case1,
      label='case1.wav after filtering with firwin')

w, h = freqz(firwin_case1, worN=8000)
plot_freq_response(w,
                   h,
                   fs_case1,
                   firwin_cutoff_case1,
                   title="Frequency response of the firwin low-pass filter")

plot_spectrogram(firwin_sound_case1,
                 fs_case1,
                 title="Spectrogram of case1.wav after filtering with firwin")

# Design an FIR filter using the formulas in my slides (i.e. do not use scipy.signal)
FIR_cutoff_case1 = .35  # I arrived at this cutoff frequency by trial and error
FIR_N = 200

FIR_case1 = FIR_cutoff_case1 * np.sinc(
    FIR_cutoff_case1 * np.arange(-FIR_N, FIR_N))
FIR_windowed_case1 = np.multiply(FIR_case1, get_window('hann', 400))
FIR_sound_case1 = convolve(FIR_windowed_case1, input_sound_case1)

# Plot the resulting spectrogram, frequency response, and play the resulting sound
sound(FIR_sound_case1,
      rate=fs_case1,
      label='case1.wav after filtering with FIR')

w, h = freqz(FIR_windowed_case1, worN=8000)
plot_freq_response(w,
                   h,
                   fs_case1,
                   FIR_cutoff_case1,
                   title="Frequency response of the FIR lowpass filter")

plot_spectrogram(FIR_sound_case1,
                 fs_case1,
                 title="Spectrogram of case1.wav after filtering with FIR")

# Design an IIR Butterworth filter using scipy.signal's routines
butterworth_cutoff_case1 = 5000  # I arrived at this cutoff frequency & order by trial and error
butterworth_order_case1 = 12
b, a = butter(butterworth_order_case1,
              butterworth_cutoff_case1 / (fs_case1 / 2))

# Filter the input sound
butterworth_sound_case1 = lfilter(b, a, input_sound_case1)

# Plot the resulting spectrogram, frequency response, and play the resulting sound
sound(butterworth_sound_case1,
      rate=fs_case1,
      label='case1.wav after filtering with Butterworth')

w, h = freqz(b, a)
plot_freq_response(w,
                   h,
                   fs_case1,
                   butterworth_cutoff_case1,
                   title="Frequency response of the IIR Butterworth filter")

plot_spectrogram(
    butterworth_sound_case1,
    fs_case1,
    title="Spectrogram of case1.wav after filtering with Butterworth")
case1.wav Your browser does not support the audio element.
case1.wav after filtering with firwin Your browser does not support the audio element.
case1.wav after filtering with FIR Your browser does not support the audio element.
case1.wav after filtering with Butterworth Your browser does not support the audio element.
In [ ]:
# Load the sound, play it, and plot it's spectrogram
fs_case2, input_sound_case2 = wavfile.read("./data/case2.wav")
sound(input_sound_case2, rate=fs_case2, label='case2.wav')

plot_spectrogram(input_sound_case2, fs_case2, title="Spectrogram of case2.wav")

# Design a filter that fixes things

# Show me the resulting spectrogram and sound
case2.wav Your browser does not support the audio element.
In [ ]:
# Load the sound, play it, and plot it's spectrogram
fs_case3, input_sound_case3 = wavfile.read("./data/case3.wav")
sound(input_sound_case3, rate=fs_case3, label='case3.wav')

plot_spectrogram(input_sound_case3, fs_case3, title="Spectrogram of case3.wav")

# Design a filter that fixes things

# Show me the resulting spectrogram and sound
case3.wav Your browser does not support the audio element.
In [ ]:
# Load the sound, play it, and plot it's spectrogram
fs_case4, input_sound_case4 = wavfile.read("./data/case4.wav")
sound(input_sound_case4, rate=fs_case4, label='case4.wav')

plot_spectrogram(input_sound_case4, fs_case4, title="Spectrogram of case4.wav")

# Design a filter that fixes things

# Show me the resulting spectrogram and sound
case4.wav Your browser does not support the audio element.

Part 2. Designing a simple equalizer¶

For this part we will design a simple graphic equalizer. We will do so using a more straightforward approach as opposed to a bank of filters as discussed in class.

We want to make an equalizer which contains six bands with center frequencies at 100Hz, 200Hz, 400Hz, 800Hz, 1600Hz and 3200Hz. Your equalizer function will take two inputs, one for the input sound and a 6-element gain vector that will indicate how much to boost or suppress each frequency band. Use the scipy.signal.firwin2 function to design a filter that has the desired characteristics. For various settings of the gain vector, use the scipy.signal.freqz command to plot the response of the filter and verify that it behaves as indicated. Experiment with various filter lengths and see which works best.

Once you figure that out, design a graphic equalizer with as many bands as you like (and arbitrary center frequencies as well), and use it to solve the problems in part 1 again. The only thing that should be different in the EQ for each recording should be the gains for each band. Play the output sounds, and show the spectrograms, see how they compare with your previous solutions.

Optional extra credit (+1pt): Use ipywidgets to make interactive sliders and process an audio stream and play it from the speakers in real-time (either from mic input, or just stream audio from disk).

In [ ]:
# Design an equalizer function
def equalizer(input_sound, gains):
    # YOUR CODE HERE
    raise NotImplementedError()


# Show its response with various gain settings
# YOUR CODE HERE
raise NotImplementedError()

# Show how it can denoise the examples in part 1
# YOUR CODE HERE
raise NotImplementedError()

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In[7], line 9
      4     raise NotImplementedError()
      7 # Show its response with various gain settings
      8 # YOUR CODE HERE
----> 9 raise NotImplementedError()
     11 # Show how it can denoise the examples in part 1
     12 # YOUR CODE HERE
     13 raise NotImplementedError()

NotImplementedError: